%--------------------------------------------------------------------------
% 脚本名称：pipeline2_3DPSF
% 用途：采用拍摄的清晰噪声图校正得到潜在清晰图像，在某一个聚焦波段位置上重建PSF
% 输出：在输出文件夹中输出PSF图像，如果有151个波段，则有151张图像。且有一个
%       mat文件保存各波段的模糊核
% 作者：刘昀卓 522543417@qq.com
% 时间：2022/1/24
%--------------------------------------------------------------------------

clc;clear;close;
% 添加必要库
addpath(genpath('../utils'));

%-----------------------------指定文件路径及参数---------------------------
SAVE_FOLDER = '结果'; % 结果保存路径
INPUT_FOLDER = 'F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据'; % 输入文件夹路径
FOCUSSHIFT_WAVELENGTH = 550; %nm 重建的三维PSF是哪个波长的聚焦位置
FOCUS_RADIUS = 30; %像素 聚焦波段圆点图像半径
DEFOCUS_WAVELENGTH = 650; %nm 模糊核大小估计所用的离焦波段
DEFOCUS_RADIUS = 48; %像素 离焦波段圆点图像半径
MAX_PSF_SIZE = 121; % 最大的模糊核估计尺寸 也是填充显示的PSF图像的尺寸
WAVELENGTH_RANGE = 420:10:720; % 要重建的波长范围
EXCEL_PATH = 'F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\版画采光谱图像位置及曝光时间数据.xlsx';% 输入的各波段聚焦位置EXCEL表
% RECTANGLE_ROI = [670,200,720,540]; % 用于裁剪图像的矩形框
savearray_path = ""; % 为空代表重新从新开始重建
cb_rows = 15;
cb_cols = 20;
rect = [278,187,174,125];
%--------------------------------------------------------------------------

%% 1.圆点目标物拍摄图测量其半径(相似原理确定核R大小)
% 1.1 聚焦搜索半径
focus_radius = FOCUS_RADIUS;
% 1.2 离焦搜索半径
defocus_radius = DEFOCUS_RADIUS;
% 1.3 通过相似原理得到比例放大系数 来为后续核估计大小
ExcelFiledata = xlsread(EXCEL_PATH,'Sheet1');% 从Excel表中导出所有数值数据
All_Wavelength = ExcelFiledata(1:end,1)';% 取出All_Wavelength,拍摄波段
All_Position = ExcelFiledata(1:end,2)';% 取出All_Position,聚焦位置
focus_shift =  All_Position(All_Wavelength==550); % 550nm聚焦位置
defocus_shift =  All_Position(All_Wavelength==650); % 650nm离焦位置
magnify = abs(focus_radius-defocus_radius)/abs(focus_shift-defocus_shift); 

%% 2.得到聚焦位置的cb图和角点
posI = im2double(imread(fullfile(INPUT_FOLDER, 'cb' ,'光谱图像',[num2str(FOCUSSHIFT_WAVELENGTH),'nm(',num2str(All_Position(All_Wavelength==FOCUSSHIFT_WAVELENGTH)),'mm).tif'])));
[fixedPoints, ~, ~] = detectCheckerboardPoints(posI);

%% 3.逐波段重建PSF
savepath = SAVE_FOLDER;
imgpath = fullfile(savepath,'PSF图片');
img1path = fullfile(savepath,'去模糊图片');
if ~(exist(savepath,'dir')==7)
   mkdir(savepath);
end
if ~(exist(imgpath,'dir')==7)
   mkdir(imgpath);
end
if ~(exist(img1path,'dir')==7)
   mkdir(img1path);
end

psf_size = MAX_PSF_SIZE;
% 加载保存的savearray
if savearray_path == ""
    savearray = cell(length(WAVELENGTH_RANGE),1);
else
    load(savearray_path)
end

for i=1:length(WAVELENGTH_RANGE)
    
    % 4.1 读入未裁剪b_img与未裁剪未放大率校正u_img
    wavelength = WAVELENGTH_RANGE(i);
    img_noise = im2double(imread(fullfile(INPUT_FOLDER, 'noise1' ,[num2str(FOCUSSHIFT_WAVELENGTH),'nm聚集位置所有中心波长图像'],[num2str(wavelength),'nm(',num2str(All_Position(All_Wavelength==FOCUSSHIFT_WAVELENGTH)),'mm).tif'])));
    img_noise2 = im2double(imread(fullfile(INPUT_FOLDER, 'noise1' ,[num2str(wavelength),'nm聚集位置所有中心波长图像'],[num2str(wavelength),'nm(',num2str(All_Position(All_Wavelength==wavelength)),'mm).tif'])));
    
    % 4.2 该波段聚焦位置放大率校正得到u_img
    LCTFI= im2double(imread(fullfile(INPUT_FOLDER, 'cb' ,'光谱图像',[num2str(wavelength),'nm(',num2str(All_Position(All_Wavelength==wavelength)),'mm).tif'])));
    [movingPoints, ~, ~] = detectCheckerboardPoints(LCTFI);
    tform = fitgeotrans(movingPoints,fixedPoints,'nonreflectivesimilarity');
    u_img = imwarp(img_noise2,tform,'OutputView',imref2d(size(img_noise2)));  
    
    % 4.3 进行PSF估计
    this_shift = All_Position(i);
    R = magnify*abs(this_shift-focus_shift); % 模糊半径
    k_size = ceil(3*R);
    
    % 将k_size变成奇数
    if mod(k_size,2)==1
        k_size = k_size+2;
    else
        k_size = k_size+1;
    end
    k_size = max(min(psf_size,k_size),15);
    
    % 裁剪
    b_img = img_noise; 
    u_img = u_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);
    b_img = b_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);

    k_init = zeros(k_size); k_init(ceil(k_size/2), ceil(k_size/2)) = 1;   

    tic
    psf = PSF_estimation(k_init,u_img,b_img,R);
    toc
    % 在填充前保存psf
    savearray{i} = psf;
    
    % 4.4 对PSF进行边缘零填充
    psf = padarray(psf,[(psf_size-k_size)/2,(psf_size-k_size)/2],'both');
    figure('units','normalized','outerposition',[0 0 1 1])
    imshow(psf,[]);
    filename = fullfile(imgpath,[num2str(wavelength),'nm.png']);
    exportgraphics(gcf,filename);
    pause(1);
    close;
    
    % 4.5 保存数组(每隔10次)
    if mod(i,10)==0
        save(fullfile(savepath,'PSF_550nm.mat'),'savearray');
    end
end
